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ABSTRACT 


This  report  summarizes  a  study  of  the  use  of  homomorphic  filtering 
and  cepstrum  analysis  for  processing  of  communication  signals  which  have 
been  corrupted  by  multipath  effects.  The  primary  objective  consisted  of 
the  development  of  computer  programs  by  means  of  which  simulated  multipath 
signals,  or  actual  received  signals  after  conversion  to  digital  form  -  can 
be  processed  so  as  to  reduce  or  eliminate  the  distorcing  effect  of  multi- 
path.  This  objective  has  been  met  for  simple  discrete  multipath.  Exten¬ 
sions  to  more  complicated  multipath  structures  are  possible.  Test  results 
have  been  obtained  using  synthesized  signal  sequences  that  simulate  actual 
communication  signals  of  a  particular  type,  namely,  16-char.nel  PSK  signals 
as  they  would  be  generated  by  an  HF  data  modem.  These  test  results  dem¬ 
onstrate  that  multipath  characteristics  can  be  identified  with  reasonable 
accuracy,  for  a  wide  range  of  multipath  parameter  values,  and  multipath 
distortion  can  therefore  be  reduced  signif  ica..tly.  Possible  practical 
implementations  are  discussed.  These  deviate  from  homomorphic  filtering 
but  utilize  the  cepstrum  as  a  means  for  identifying  the  multipath 


structure. 
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1.  INTRODUCTION 


This  report  summarizes  a  study  of  the  use  of  homomorphic  filtering  and 
cepstrum  analysis  as  described  by  Oppenheim,  Schafer  and  Stockham  [1]  for 
processing  of  communication  signals  which  have  been  corrupted  by  multipath 
effects.  The  primary  objective  consisted  of  the  development  of  computer 
programs  by  means  of  which  simulated  multipath  signals,  or  actual  received 
signals  after  conversion  to  digital  form  ~  can  be  processed  so  as  to  reduce 
or  eliminate  the  distorting  effect  of  multipath.  This  objective  has  been 
met  for  simple  discrete  multipath.  Extensions  to  more  complicated  multipath 
structures  are  possible.  Test  results  have  been  obtained  using  synthesized 
signal  sequences  that  simulate  actual  communication  signals  of  a  particular 
type,  namely,  16-channel  PSK  signals  as  they  would  be  generated  by  an  HF  data 
modem.  These  test  results  demonstrate  that  multipath  characteristics  can  be 
identified  with  reasonable  accuracy,  for  a  wide  range  of  multipath  parameter 
values,  and  multipath  distortion  can  therefore  be  reduced  significantly. 
Possible  practical  implementations  are  discussed.  These  deviate  from  homo¬ 
morphic  filtering  but  utilize  the  cepstrum  as  a  means  for  identifying  the 
multipath  structure. 

This  report  is  arranged  into  major  Sections  as  follows: 

Section  2  -  Homomorphic  Filtering  and  the  Cepstrum.  The  basic  non-linear 
signal  processing  ideas  which  lead  to  the  so-called  "cepstrum"  are  briefly 
surveyed,  in  terms  of  continuous  signal  representation.  It  is  shown  that  if 
two  signals  are  convolved,  their  cepstra  are  summed.  Multipath  is  therefore 
an  additive  phenomenon  in  the  cepstrum. 

Section  j  -  Discrete.  Finite-Duration  Cepstral  Analysis.  The  discussion  in 
Section  2  in  terms  of  continuous  signals  is  here  extended  to  discrete 
(sampled)  signals  which  are  observed  over  a  finite  observation  interval. 
Several  deviations  from  the  continuous  analysis  are  pointed  out. 

Section  4  -  The  PSK  Signal.  In  this  Section,  the  16-channel  PSK  signal  is 
described  which  was  simulated  for  used  in  all  performance  tests. 

Section  5  -  The  Computer  Program  The  program  is  described  which  has  been 
evolved  to  investigate,  detect  and  remove  multipath  effects  from  simulated 
multipath  signals.  Sufficient  detail  is  presented  so  that  the  various  pro¬ 
gram  functions  and  options  can  be  understood  and  the  program  can  be  run. 
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Section  6  -  System  Performance.  The  effect  of  errors  in  the  estimates  of 
multipath  parameters  on  the  removability  of  the  multipath  distortion  is 
analyzed.  On  the  basis  of  this  analysis,  estimates  of  system  performance 
for  simple  discrete  multipath  are  obtained,  neglecting  additive  noise.  It 
appears  that  with  good  reliability,  multipath  interference  can  be  reduced 
to  a  20  db  signal- to- interference  ratio  over  a  wide  range  of  multipath 
amplitude  and  delay. 

Section  7  -  Discussion.  Possible  improvements  of  the  present  system  are 
suggested,  and  possible  forms  of  practical  implementation  are  considered. 
More  experimental  data  is  needed  in  order  to  arrive  at  adequate  specifica¬ 
tions  for  a  practical  system. 


2 .  HOMOMORPHIC  FILTERING  AND  THE  CEPSTRUM 


2.1  The  Logarithmic  Power  Spectrum 

The  logarithmic  spectrum  was  described  by  Bcgert,  Healy  and  Tukey  [2]  as 
an  approach  to  the  separation  of  echoes  from  a  received  signal.  Thus, 
consider  a  received  sigr-1  of  the  form 

y(t)  =  x ( t )  +  ax(t-x) ,  \a\  <_  1,  (2-1) 

where  x(t)  is  the  transmitted  signal.  Then  it  may  be  desired  to  identify  the 
parameters  of  the  echo,  that  is,  a  and  r;  or  it  may  be  desired  to  extract 
an  undistorted  version  of  x(t). 

The  autocorrelation  function  of  the  (real)  signal  x(t)  is 

A 


ft  00 

(ft^(£)  a  |  x(t)  x(t+£)dt,  finite  energy  case 


> 


x 

lim  -Tfz  I  x(t)  x(t+£)dt,  infinite  energy  case  (2-2) 

T*oo  -1  J_x  J 


If  (ft^ (O  has  a  narrow  peak  and  is  small  elsewhere,  autocorrelation  analysis 

of  the  received  signal  y(t)  allows  determination  of  the  echo.  Expressing 

(ft  (O  in  terms  of  (ft  (?) , 
yy  xx 


^yy^)  =  *sv/0  +  a2  +  <*  ^xx(f.-1>  +  a^XX^+T^> 


(?’  J) 


the  magnitude  of  the  central  peak  is  seen  to  be  (1+a  )  (ft  (0)  +  2a  (ft  (t)  « 

2  xx  xx 

(1+a  )  (R  (0)  ,  and  the  magnitudes  of  the  peaks  at  ±j  are  a  (ft  (0)  +  (l+<*2i(ft  (r)+ 

XX  XX  XX 


XX 


Fig.  1  Sketch  of  Autocorrelation  Function  with  Echo  Component 

In  many  cases  of  practical  interest,  however,  autocorrelation  analysis  may 
give  inconclusive  results,  This  is  particularly  true  if  only  relatively 
short  signal  samples  are  available  for  pro'-  '_-sj.ng. 
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Bogert,  Healy  and  Tukey  therefore  examined  the  power  spectrum  of  the 
received  waveform  in  order  to  find  a  way  of  isolating  the  echo  component. 
This  led  to  methods  of  processing  the  received  signal  which  involved 
taking  the  logarithm  of  the  power  spectrum.  (See  Table  1). 


Table  1. 

Signal 

F.T. 

Power  Spectrum 

Log  Spectrum 

x(t) 

«K(f> 

§(f)  -  |^|2 

£  Cf) 

X 

^XX  X 1 

XX 

y(t) 

(l+e-^f 

>  y,<f> 

A(£)-gxx(£) 

In  A(f )  +£  (f) 

XX 

where  A(f) 

2 

*  (1+a  +2a  cos 

2trf  t) 

Thf  function  A(f)  in  Table  1,  and  its  logarithm,  have  the  form  shown  in 
Fig.  2.  In  the  log  spectrum  the  nearly  sinusoidal  ripple  of  In  A(f)  is 
superimposed  on  £  ,  the  log  spectrum  of  x(t).  The  form  of  £  depends 


on  x(t),  but  come  general  comments  can  be  made.  Since  the  power  spectrum 
generally  decays  torard  zero  with  increasing  |f|,  the  log  spectrum  goes 
toward  large  negative  values  with  increasing  |f|.  Small  ripples  in  §  in 

XX 

a  region  where  §  is  very  small  will  result  in  very  large  ripples  in  £  . 

XX  XX 

Depending  on  the  units  used  in  computing  Q  ,  it  is  possible  that  &■  (f)  <  1 

XX  XX 

for  all  f.  In  that  case,  j£^(f)  is  everywhere  negative  -  i.e.,  it  has  a 
large  negative  "d.c.  component."  The  latter  is  also  present  due  to  the 
negative-going  trend  with  large  jf|,  mentioned  before.  A  typical  § and 
corresponding  £  migh  appear  as  in  Fig.  3.  There  remains  now  the  task  of 

XX 

extracting  the  periodic  ripple  from  £  ,  which  Is  basically  a  problem  of 

detecting  a  sine-wave  in  the  presence  of  noJse.  T^s  can  be  accomplished 
e^her  through  autocorrelation  or  spectral  analysis,  both  of  which  are 
discussed  by  Bogert,  Healy  and  Tukey.  Since  the  "signal"  is  now  a  frequency 
function,  they  speak  of  its  cepstrum .  rather  than  spectrum,  and  of  its 
pseudo-autocovariance,  rather  than  autocovariance.  The  cepstrum  has  been 
included  in  the  continuation  of  Table  1.  As  can  be  seen, 


Table  1  (continued) 

Log  Spectrum 

Cepstrum  (Power  spectrum  of  log  spectrum) 

£  (f) 

XX 

£  (f)  +  2t.  cos  Zirf t 

XX 

( 9  (t’) 

XX 

Q  ( t ’ )  +  a^6(t’ -x)  +  n^5(t'+i) 

XX 
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the  delay  component  in  y(t)  has  the  effect  of  superimposing  on  the  cepstrum 
of  x(t)  a  pair  of  impulses, located  at  t'  =  t  t  (see  Fig  4). 


Fig.  4  Sketch  of  a  typical  signal  cepstrum  with  superimposed  impulse 
pair  due  to  delay  component. 

Due  to  the  fact  that  lnA(f)  is  not  purely  sinusoidal,  there  must  alto  be 
small  contributions  at  multiples  of  t 

2 . 2  Hoiik  aorphic  Deconvolution 

A  different  approach  to  the  identification  and  removal  of  echoes  was 
pursued  by  Oppenhei.m,  Schafer  and  Stockham  [1],  based  on  the  theory  of 
homomorphic  filtering  developed  by  Oppenheim  [3].  This  calls  for  a  non¬ 
linear  processing  system  of  the  form  shown  in  Fig.  5. 


Fig.  5  Block  diagram  of  Nonlinear  Processing  System  for  Homomorphic 
Filtering. 

The  nonlinear  blocks  D  and  D  ^  perform  homomorphic  transformations  of  a 

signal  space.  In  the  application  at  hand,  this  means  that  D  transforms  a 

“convolutional  space"  into  a  "sum  space";  that  is,  two  signals  which  are 

convolved  at.  the  input  result  in  images  that  are  added.  Thus,  if  two 

different  signals  x  x^  applies  to  D  produce  outputs  x^>  respectively, 

then  the  convolution  of  x^  and  xa2>  denoted  xa^*xa2’  Pr°duces  the  output 

x,  ,-tx,  „  (Fig.  6).  This  particular  transformation  is  also  called  homomorphic 
bl  b2 

deconvolution. 

In  order  to  accomplish  homomorphic  deconvolution,  D  must  perform  three 
operations  in  succession:  Fourier  transformation  of  the  input  signal ,  complex 
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logarithm  of  the  resulting  transform,  followed  by  inverse  Fourier  trans¬ 
formation. 


al 


a2 


- >. 

Sl  +  S2 


Fig.  6  Deconvolution 

D*  is  the  inverse  of  D.  In  the  configuration  of  Fig.  5,  it  may  then  be 
possible  to  deconvolve  two  input  signals,  through  linear  filtering  (block 
L  in  Fig.  5)  and  then  return  to  the  input  signal  space  via  the  D 
transformation 

When  applied  to  convolved  signals,  homomorphic  filtering  becomes 
very  similar  to  the  generation  of  the  cepstrum.  In  fact,  the  output  obtained 
after  the  homomorphic  deconvolution  operation  (block  D  in  Fig.  5)  may  be 
called  the  complex  cepstrum  of  the  input,  since  the  complex  logarithm  is 
involved.  The  details  are  shown  in  Table  2. 


Table  2 

Signal 

F.T. 

Complex  Logarithm 

Complex  Cepstrum 

x(t) 

Tx«> 

£  (f)  =  cln  <T  (f) 

X  X 

<9  (f)  =  inverse  F.T. 
X  of  £x(f) 

y(t) 

Tx(f) 

cln(.'+ae  ^^  ‘  )+.£  (f ) 

see  below 

First  it  is  clear  that  the  logarithmic  power  spectrum  of  the  preceding 
Section  has  a  simple  relation  to  the  complex  logarithm  of  the  Fourier 
Transform: 

£  (f)  =2  He  (f)  (2-4) 

XX  X 

The  cepstrum  of  Bogert,  Healy  and  Tukey  is  therefore  twice  the  even  part  of 
the  complex  cepstrum: 

Q  (f)  =  2  Ev  <2  (f)  (2-5) 

XX  X 
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(For  an  even  function,  forward  and  inverse  Fourier  transforms  produce  the 
result).  Next  the  factor  (1+ae  ^1T^T)  in  ^(f)  has  to  be  considered. 


same 


The  locus  of  this  function  in  the  complex  plane  is  a  circle  centered  at 
1  +  jO,  with  radius  a,  so  that  cln(l+ae  ^  71  l)has  the  form  shown  in  Fig.  7. 


(b) 


Fig,  7  Th .  Function  1+ae  ^  17-1 . 

(a)  Locus  in  complex  plane 

(b)  Waveforms  of  real  and  imaginary  parts. 

These  approximately  sinusoidal  ripples  are  the  effect  of  the  echo  and 
are  superimposed  on  the  real  and  imaginary  parts  of  The  complex 

cepstrura  of  y(t)  is  therefore  of  the  form 


Q  {t')  “  v?  (t*)  +  j^.T.  of  cln(l+-e"j2lTfT)  . 


(2-6) 


- i 2v£t 

Since  cln(l+-e  J  )  is  a  periodic  waveform,  its  inverse  transform  consists 


of  spikes  at  the  fundamental  and  harmonics  (on  tire  t'-axis).  Coefficients 
of  the  impulses  appealing  in  €'  (t1)  are  obtained  by  expanding  the  logarithm 


cln(l+ae  ^'^T) 


ae 


-j4nfT  +  a_  e“ j 6tt f t 


(2-7) 
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The  inverse  transform  then  results  in 


2  3 

a6(t'-T)  -  6(t’-2t)  +aj6(t'-3i)  -  ....  (2-8) 


Representative  coefficient  values  are  shown  in  Table  3. 

Table  3 


Coefficients 

of  (2-8) 

.1 

.1  , 

.005, 

.00033 . 

.5 

•  5  , 

.125, 

.042,  _ 

.9 

.9  , 

.405, 

.243, 

.99 

.99, 

.490, 

.323,  - 

Since  the  signal  y(t)  in  (2-1)  is  obtained  by  passing  x(t)  through  a 
channel  with  impulse  response  h(t)  =  6(t)  +  a6(t-i),  i.e., 

y(t)  *  x(t)  *  1 6 ( t)  +  a6(t-r)], 

the  expression  (2-8)  represents  C.(t'),  the  complex  cepstrum  of  h(t). 

a  comb  filter  might  now  be  used  to  eliminate  (?  ( t * )  from  Q  (t'), 

n  y 

leaving  a  relatively  undistorted  version  of  (?  (t').  Application  of  the 

-1  X 
D  operation  should  then  yield  the  reconstructed  signal  x(t). 

2 . 3  More  Complicated  Multipath  Structure 

If  the  received  signal  y(t)  is  of  a  different  form  than  has  been 
assumed  so  far,  the  situation  is  considerably  more  complicated. 

In  the  case  of  two  delay  components,  the  following  expressions  apply: 

y (t)  =  x(t)  +  ax(t-x^)  +  6x(t-T2)  ,  (2-9) 

tF  (f)  =  (1  +  arj27tfTl  +  6e_j2"fT2)  ®F(f)  (2-xO) 

y  x 

£(t)  =  cln(l4de"j27lfT1  +  6e"j2nfT2)  +  £  (f)  (2-11) 

y  x 

Here  the  nature  of  the  complex  cepstrum  can  be  inferred  without  resorting 
to  detailed  calculations.  If  the  effect  of  the  complex  logarithm  in  (2-11) 
is  neglected  and  merely  the  inverse  transform  of  the  term 
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(14ue  ^1T^T1  +  3e  ^1T^T2)  is  taken,  this  would  introduce' in  £?  (t1)  the 
impulses  6 ( t * )  +  a6(t'-t^)  +  g6(t-T?).  Due  to  the  non-linearity  of  the 
logarithm,  however,  harmonics  and  intermodulation  terms  are  produced. 
Thus,  the  second-order  contributions  appear  at  2t^,  and  +  t  . 

More  specifically  (2-7)  is  now  replaced  by: 

— j  27tf  t  -j2TTfT_  -J2ttxt  -j2nfr? 

cln(l+ae  +  Be  )  =  (cie  +  Be  ) 


It) 


The  term 


%«> 

—  in  effect  represents  a  random  phase  and  amplitiv 

modulation  of  ae  ^7r^T .  This  modulation  appears  in  the  cepstrum  domain 
in  the  form  of  noise  sidebands  about  the  spike  pattern  which  identifies 
the  multipath  characteristics.  The  spike  pattern  itself  is  not  affected, 
however,  as  long  as  inf ini te- time,  continuous  analysis  is  used.  Acceptable 
behavior  of  ^(f)  is  also  assumed. 
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3 .  DISCRFTE.  FINITE  DURATION  CEPSTRAL  ANALYSIS 


Whan  working  with  finite  duratijn  sample  sequences,  as  is  the  case 
when  using  'Fast  Fourier  Transform"  (FFT)  techniques,  the  analysis  of  the 
previous  chapter  is  not  directly  applicable.  Some  of  the  differences  will 
be  explained  in  this  chapter. 

3.1  The  Discrete  Cepstrum 


In  Table  4,  a  comparison  between  the  continuous  and  discrete  finite- 
duration  cepstral  analysis  is  drawn  up.  (Braces  denote  sequences  of 
numbers. ) 


Table  4.  Continuous 

and  Discrete  Cepstral  Analysis 

continuous 

discrete 

The  signal 

X(t)  i  -oo  <  C  <  oo 

{x^}  k  =  0, . . . ,N-1 

where  the  k1"*1  sample  value  x^xC^-) 

The  spectrum 

T(f) 

,  ,  ( i  N-l  -12itnk  i 

‘is  |0v  "  |  • 

n  =  0, . . . ,N-1. 

The  log 

spectrum 

{ l  (n)  t  =  cln{  d  (n)  } 

X  X 

The  complex 
cepstrum 

<?xU') 

.  2rrnk 

[n-1  J  N  | 

(c  (k)}®|^2^.  (n)e  |k=0,...,N-l 

X  ln-0  X  j 

The  discrete  spectrum  {dx(n)>  in  general  is  not  ju.->t  a  sampler  version  of 
CF  (f)»  but  includes  an  aliasing  effect,  (if  x(t)  is  not  bandlimited  to 
frequencies  below  ~)and  the  effect  of  truncation  in  the  time  domain. 

The  precise  relation  between  sample  values  of  and  the  sequence  {£  }  is 


12 


somewhat  obscure,  because  the  effects  of  aliasing  and  time-domain  trunca¬ 
tion.  which  are  present  in  {d  },  have  become  intermixed  with  the  sample 
values  of  <j£  .  The  sa'ie  is  true  for  the  relation  between  <2  and  {c  }. 

X  XX 

However,  these  effects  are  small  if  a  large  record  is  processed,  and  if 
sampling  is  sufficiently  rapid  so  that  little  aliasing  arises.  The  dis¬ 
crepancy  due  to  truncation  in  time  disappears  if  the  signal  x  is  periodic, 
with  period  T. 

Discrete  cepstral  analysis  can  also  be  discussed  in  terms  of  z- 
transforms.  The  z-transform  representation  of  the  signal  sequence  {-x^} 
is  defined  as 

N-l 

Z  (z)  *=  X  v"*’  (3~1) 

x  k-0 

so  that  j2nn 

(dx(n)}  -  |  {  Z;c(e  N  )}.  (3-2) 

This  is  a  sequence  of  N  values  of  3^. (z)  which  occur  at  equally  spaced 
points  along  the  unit  circle  in  the  z-plane,  multiplied  by  1/N.  Then 

j  2nn 

{£  (n)}  *  cln{d  (n) }  »  -InN  +  { cln  Z  (e  ^  ))  (3-3) 

XX  x 

But  cln  Zx  can  itself  be  regarded  as  the  z-transform  of  some  sequence. 
Thus,  let 

Z  (z)  “  -NlnN  +  N  cln  Z  (z) ,  (3-4) 

C  X 

X  ;2wn 

so  that  it  (n)}  *  Z  (e  N  )}.  Comparison  with  (3-2)  shows  that 
x  N  c 

(l  (n) }  is  then  the  DFT  of  the  sequence  of  numbers  whose  z-transform  is 

x  _k 

2  (z)  c  (k)z  .  In  other  words, 

Cx  k*0  X 

(1  (n)}  is  the  DFT  of  the  complex  cepstrum  (c  (k) ) , 

X  x 


",  * 


3.2  The  Multipath  Cepstrum  -  Simplest  Case 


In  order  to  see  how  the  effect  of  multipath  appears  in  the  discrete 
cepstrum,  we  consider  first  the  case  of  a  periodic  signal  with  period  T, 
and  a  multipath  delay  which  equals  an  integral  number  of  sampling  intervals. 
Thus,  v;e  have  a  transmitted  signal  x(t)  =  x(t-T)  and  a  received  signal 
y(t)  as  in 

y(t)  =  x { t)  +  a(t-x) , 


but  with 


T 


rT 

N’ 


where  r  is  any  integer  from  1  to  N-l. 


(3-5) 


Since  x  is  periodic,  <3r  is  discrete.  If,  in  addition,  x  is 

X  N 

strictly  bandlimited  to  frequencies  less  than  — ,  then  the  coefficients 

of  <Jr  within  the  band  limits  coincide  with  (d  (n' (see  [10].) 

X  X 

The  assumption  that  x  is  strictly  bandlimited  is  not  needed  here, 

however,  and  it  will  be  understood  that  the  discretf  spectrum  {  d^(n) }may 

differ  from  . 

x 

In  place  of  the  F.T.  of  y  (Table  2),  the  DFT  of  y  is  now  considered: 

{d  (n)>  -  (d  (n)(l  +  ue  J  N  )},  n=0,l,..,S-l  (3-6) 

y  ^ 


nr 

-  j  2ir—“ 

{  £  (n)  }  =  {£  (n) }  +  {cln(l  +  ote  )} 

y  x 


nr 


-1  'J  2'*N 

(c  (k)}  =>  {c  (k)l  +  DFT  (cln[l  +  ae  ]} 

y  ^ 


(3-7) 

(3-8) 


The  second  term  in  (3-8)  gives,  using  (2-7): 


,  .  „  nr  2Ttnk  -j 

,  -  i  -  v*  J  IT" 

{  cm(k)}  c3n(^-  +  ae  )  e 

(n=0 


i 


■ft  (-£** 

[n-0  '  p=l  1 


,  ~  nr .  ,  2^nk 

-a)P  ~12rpin  j  " 
— e 
n 


e 


N 


1 
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rK-a).1 


for  k  =  pr  (p  =  1,2,...)  (mod  N).| 


(3-9) 


0  otherwise 


Alternately,  the  sampled  representation  of  v(t)  can  be.  expressed 
in  terms  of  its  z- transform: 


Z  (z)  *  Zx(z)  +  otzr  2x(z)  *  2^(z)(l  +azr). 


(3-10) 


Then 


2  (z)  -  -NlnN  +  N  cln  2  (z) 

Cy  '' 


*  -NlnN  +  N  cln  2  (z)  +  N  cln  (l+azr) 


-  2 :  <«)  -  s  Z. 


is)  ~rP 


(3-11) 


p=l 


The  last  term  is  the  z-transform  of  the  multipath  cepstrum  {c  (k)}, 

m 

so  that  (3-9)  is  again  obtained.  In  other  words,  superimposed  on  the 
discrete  cepstrum {c^j appears  a  multipath  pattern  such  as  shown  in  Fig.  8 
for  the  particular  case  a  =  0, j,  r  =  5,  N  =  512 


3CC  4 


20C 


lOc 


0  4 


-  ICO 


■o  o  o  o  L~ e-e-a-a— r-  «-o--e-o - o-o  t>  o  -  ■>  A 

)  2  3  4  i  t,  /  8  9  13  il  12  I)  M  !*-  !6  17 

O 


Fig.  8  Typical  Discrete  Multipath  Cepstrum 
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3.3  Multipath  Cepstrum  -  Delay  is  not  an  integral  number  of  sampling 
intervals 

Next,  an  arbitrary  multipath  delay  t,  is  considered  but  x  still 

periolic  with  period  T«  Here  x(t)  will  at  first  be  assumed  bandlimitad 

N 

to  frequencies  below  —  .  The  received  signal  is 

y(t)  =  x(t)  +  ux(t-t) , 

or  in  terms  of  sample  values, 

{yfc}  -{xk  +  ax(T)k}  (3-12) 


where  the  x,  are  samples  of  a  delayed  version  of  x  (delayed  by  t  sec.) 
v  T  )  K 

Since  x  is  bandlimited,  the  x^k  can  be  expressed  exactly  in  terms  of 
the  xk's,  using  the  sampling  expansion  of  x(t): 


x(t)  =  E  xk  sine  (“■  -  k) 
k®-® 


(where  x.  ... 

k+N 


x^ ,k  —  • • . —1 , 0 , 1 , . • . ) • 


Thus, 

x (t ) k  =  X^Nk  “  T)  =  E  x(p^)  sine  (k  -  t  ~  -  p) .  (3-13) 

'  '  p=-oo 

Since  x  is  taken  to  be  periodic  with  period  T,  only  N  distinct  coeffi¬ 
cients  are  involved  in  (3-13) ,  which  can  therefore  be  written  in  the 
following  way  in  order  to  exhibit  that  fact: 


“(T)k 


N-l 

E 

r=— oo  p=o 


-I 


sine  (k  -  t  —  -  p  -  rN) . 


(3-14) 


Applicatic  n  of  the  TIT  to  x^k  gives 

N-l 

•  : 

P 


32  ^ 


1  4—  r- —  M  J  N 

dv  (n)  *  :  I  E  x  tine  (k  -  t  -  -  p)e 
X(T)  N  £:»  p  T 


N-l  ®  N-l  M  -j  2tt^- 

*  N  E  E  xp  sinc  -  T  t  '  P  "  rN)e 

k=0  r— »  p=0  P 

(3-15) 
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The  DFT  of  {y^}  i'  f».  ref  ore 


dJr>  “  N  ^  [*k  +  Xp  Sinc(k  "  T  T  “  p)]  e  ^ 


»  N-l 


■dx<n)+»2I  z  X  sinr-x  -T  f  -  p) 


-j  2?r^- 

3  V 


(3-lfc) 


~  "  p=-co  k=0  y 

The  second  term  becomes,  upon  replacing  p  by  k  +  p: 


fZ  Z  xk+  PlncU  |  +  p) 


.  ^  nk 

-J2v 


p=-“  k=0 


°°  j  2ir^ 

«d  (n)  sine (t  ^  +  p)e  N  ; 


so  that  (3-lfc)  becomes 


d  (n)  *  d 

V 


*  00 

i  (n)  1  +  ot  y 

D=—  cc 


Oir*2,  ' 
N  J  N 

sinc(T  -  +  p)e 


(3-17) 


This  should  be  compared  with  (3-6)  ,  to  which  it  reduces  if  t  is  a 
T 

multiple  of  The  discrete  cepstrnm  then  becomes 

[  00  j  2ir^  ) 

c  (k)  *  c^tk)  +  D.F.T.  'Vein  [1  +  sine  (t  ~  +  p)e  N 

7  X  l  P=-°°  J 

-  cx(k)  +  D.F.T.-1  (aL—)-  2  (“!••••)  +3^S*-')- 


(3-18) 


The  first-order  multipath  term  is,  for  k  =  0,1,..., N-l, 

b  -l  00  M  j  2tt““— 

cml(k)  =  a  sine  (t  —  +•  p)e  e 

l!  »'.=0  p=-» 
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oo  N-l  j2rr'n'0^tH) 

X,  sine  (t  -  +  p)e 
p*-«  n=0 

Only  those  terms  for  which  p  =  -k  *  mN  are  terms  which  contribute  to  the 
result,  giving 

CO 

c  .  (k)  *  aN  y  sine  (t  *r  -  k  +  mN)  =  a  sin  tt(t  7  -  k)  cot  u(~  -  7) 
ml  — •  i  T  t  N 


aN  sine  (t  ^  -  k)  ,  for  ^  -  ^<<  1. 


(3-19) 


The  >econd-order  multipath  term  is,  for  k  =  0,1,..., N-l 

N-l  2  -  co  , 27TJEliJG!±ajL  i2rnk 

cm2(k)  =  X  %  X  X  slnc  (t  f  +  p)  sinc  (T  ^  q)e  N  e  N 

n*0  p=-oo  q=-oo 


2  00  00 

-f-  I  I 


N-2  j2ir 


p-—  ^Iloo  sinc  -  +  p)  sinc  (t  -  +  q)  2L  e 


Only  those  terms  contribute  for  which  p  +  q  =  -k  j)  mN.  Then 

2  00  CO 

Cm2^k^  3  ^2^  X  X.  sinc  (t  ^  -  k  -  q  +  mN)  sinc  (x  ^  +  q) . 
m*-”  q»-« 

The  inner  summation  is  a  discrete  convolution  of  two  sinc-functions  which 
yields  again  a  sinc-function: 


2  » _  „  t 

Cm2^k>  ”  ^2^  X  sfnc(2x  7  -  k  +  mN)  =  ~  sir.  it ( 2t“  -  k)  cot  n(~  -  ~) 


..  o.2N  .  N  ,  2x  k  . 

2~  sinc  (2t  7  ~  *0  ..  for  ~ - -  <■<  1 


(3-20) 
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In  a  similar  way  are  the  remaining  multipath  terms  approximated  by 
discrete  sine-functions.  (3-18)  therefore  becomes 


cy(k) 


(k) 


-Hi 

p»l 


(-a)P  N 

_  sinc(pT  ~  -  k)  , 


k  =  0,1,. • • »N-1 


(3- 


The  argument  of  the  sinc-function  is  to  be  understood  as  mod  N  for  each  p 


Thus  there  appears  superimposed  on  the  discrete  cepstrum  {  c^}  a 
multipath  pattern  such  as  shown  in  Fig.  9  for  the  case  ot  =  0.5, 
x  *  5.5  N  -  512: 

T'} 


Fig.  9  Discrete  Multipath 
Cepstrum;  t  not  a 
multiple  of  the 
sample  spacing. 


3.4  Multipath  Cepstrum  -  Arbitrary  Signal 

For  an  arbitrary  signal  x(t),  two  sources  of  error  enter  into 

expression  (3-20).  First,  if  x(t)  is  not  periodic  witn  period  T,  the 
T 

values  xCirk-x)  (for  k  =  0,1,...,N-1)  are  no  longer  expressible  exactly 

T 

in  terms  of  the  samples  x(k^) .  Instead  of  (3-13),  only  an  approximate 
expression  can  be  given, 

Nrl  T  N 

x(T)k  sinc(k  -  T-  -  p),  (3-22) 

T 

where  the  approximation  is  reasonable  only  for  — k  -  x  within  the  interval 

N 

(0,T),  and  not  too  close  to  the  end  points  of  this  interval.  Second, 
if  x(t)  is  not  bandlimited  to  frequencies  below  ”,  then  an  additional 
error  is  introduced  due  to  frequency-domain  aliasing . 


21) 
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The  two  effects  just  described  are  not  very  significant  if 

a)  rather  long  sample  sequences  are  processed; 

b)  the  signal  is  sampled  sufficiently  rapidly  so  that  a  negligible 

fraction  of  signal  energy  lies  outside  the  frequency  range  below 
N 

2^;  and 

c)  the  maximum  multipath  delay  t  is  much  smaller  than  T. 

Theif  requirements  have  been  assumed  to  hold  in  this  investigation.  The 
multipath  detection  procedure  described  in  Chapter  5  has  therefore  been 
based  on  Eqs.  (3-9)  and  (3-21, . 

The  discrete  analysis  extends  in  an  analogous  manner  to  more 
complicated  multipath  structurts,  such  as  describe^  in  Equ.  (2-9)  through 
(2-12)  in  terms  of  continuous  cepstral  analysis.  However,  the  multi- 
path  detection  and  cancellation  procedure  presently  programmed  and 
described  in  Chapter  5  operate  only  on  a  single  multipath  component. 

3.5  Additive  Noise 

Extension  of  Equ.  (2-16)  to  the  finite  duration  discrete  case  gives 
for  the  discrete  cepstrura  of  a  multipath  signal  with  additive  noise: 


cx(k)  + 


D.F.T.'1 


1  +  ae 


dv(n) 

dx(n) 


(3-23) 


For  convenience,  multipath  delay  is  assumed  a  multiple  of  the  sampling 
interval,  as  in  Sec.  3.2,  and  is  expressed  by  (3-5).  As  in  the 
continuous  case,  the  term  can  be  regarded  as  random  phase 


dx(n) 


and  amplitude  modulation  of 


,  which  introduces  noise  side¬ 


bands  about  the  multipath  spikes  in  the  discrete  cepstrum.  In  the 
discrete  case,  however,  because  of  the  finite  observation  time,  a  random 
error  will  also  appear  superimposed  on  the  magnitudes  of  the  multipath 
spikes.  A  detailed,  quantitative  derivation  of  thi:  effect  has  not 
been  carried  out.  However,  from  a  first-order  analysis  it  follows  that 
the  noise  does  not  bias  the  magnitudes  of  the  multipath  spikes  in  the 
real  cepstrum. 
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In  order  to  show  this,  attention  is  restricted  to  the  real 


logarithm  in  (3-23),  that  is. 


. nr 

-J2TtTT 


N 

In  1  +  ae  +  z(n) 


(3-24) 


where  z(n)  *  d^(n)/dx(n).  The  phasor  diagram  in  Fig.  10  is  helpful: 


The  phasor  with  magnitude  c  advances  with  n  in  discrete  steps  of  size  — 
For  given  n,  its  sum  with  the  unit  vector  has  been  denoted  y  in  Fig.  10, 


where 


'1  +  2a  cos  2tt— — 
N 


(3-25) 


Only  the  magnitude  of  y  is  needed  since  0,  the  argument  of  the  random  pha¬ 
sor  z(n) ,  is  assumed  uniformly  distributed  over  (0,2tt),  and  independent 
of  | z|  and  arg  y.  Aduing  to  y  the  random  phasor  z  results  in  n,  where 


(3-26) 


In  order  to  average  out  the  effect  of  0,  the  expected  value  with 
respect  to  0  can  be  computed.  This  can  be  done  as  follows,  under  the 
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assumption  that  the  noise  is  small  so  that  |~|  <  1  for  all  n: 

Y 

1  .  2 

Eeln|nj  =  E0{ In 1 y t  +  jln(l  +  2 1^|  cos  9  +  |^|  )), 

2 

where  -|ln(l  +  2|^|  cos  0  +  |— I  )  =  (Rx.  eln(l  +  |^|e^),  so  that 

1  2 

E_ln|n|  =  In | y |  +  E_{  |— !  cos  0  -  t|~|  cos  29  +  ...} 

9  tty  i-  y 

=  In | y | •  (3-27) 

The  expected  values  of  the  multipath  cepstrum  are  therefore  not 
affected  by  the  noise,  in  the  case  of  large  signal- to-noise  ratio. 

3.6  Treatment  of  the  Errors  Described  in  Sec.  3.4  as  Additive  Noise 

If  the  transmitted  signal  is  not  assumed  periodic  with  period  T, 
then  a  discrepancy  arises  due  to  the  fact  that  the  observed  portion  of 
the  delayed  signal  component,  ax(t-t),  is  not  merely  a  delayed  version 
of  the  observed  portion  of  the  transmitted  signal  component  x(t). 


Fig.  11  A  L'ignal  and  its  Delayed  Replica 
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Fig.  11  shows  a  signal  waveform  x(t),  and  its  delayed  replica, 
x ( t— t ) .  When  observation  is  confined  to  (0,T),  the  delayed  signal  as 
observed  does  not  correspond  to  the  observed  portion  of  x(t) ,  but  only 
to  x(t)  within  (0,T-x).  Also,  in  the  interval  (0,~),  the  delayed 
waveform  does  not  agree  with  the  observed  x(t)  anywhere.  This  can  be 
viewed  in  the  following  way.  Let 


x0(t) 


0  <_  t  <  T-  r 

T-t  <  t  <  T 


x1(t) 


T-t  <_  t  <  T 
0  <-  t  <  T-t 


x_2(t) 


-T  <_  t  <  0 
otherwise 


>  (3-28) 


i 
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Then  the  observed  portion  of  x(t)  can  be  expressed  in  terms  of  two 
components , 


x(t)  =  xQ(t)  +  x^t) 


(3-29) 


The  observed  portion  of  the  received  multipath  signal  can  be  written 

y(t)  -  x(t)  +  cxx(t-T)  «  xQ(t)  +  ax^Ct-T )  +  x^t)  +  ax_1(t-T).  (3-30) 

Within  the  observation  interval,  only  the  x^  component  appears  to 
have  undergone  multipath  delay.  Therefore,  when  the  received  signal 
(3-30)  undergoes  finite-duration  discrete  cepstral  analysis,  the  x^ 
and  x_^  terms  can  be  viewed  as  noise  which  does  not  carry  any  multi- 
path  information.  The  signal-to-noise  ratio  can  still  be  considered 
sm<il)  .  as  long  as  t  <<  T.  However,  the  waveform  segments  x^(t)  and 
ax_^(t-T)  arc  dependent  (in  amplitude  and  duration)  on  the  multi- 
path  characteristics.  Therefore,  their  mean  effect  on  the  multipath 
cepstrum  is  not  zero,  as  would  be  the  ~ase  for  independent  low-level 
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noise.  Instead,  the  multipath  cepstrum  appears  as  an  attenuated  version 
of  the  ideal  multipath  cepstrum,  with  an  attenuation  roughly  pro¬ 
portional  to  t/T. 

3.7  Use  of  a  Data  Window 

The  effect  just  described,  as  well  as  spectrum  spread  due  to 
truncation  of  the  time  waveform,  can  be  reduce  l  by  the  use  of  a  data 
window.  Only  one  type  of  window,  a  "Hanning"  window,  has  been  tried. 

This  consists  of  multiplying  the  first  26  data  points,  x_,  x_,..., 

1  nit  ^  ^ 

x25>  by  ^(1  -  cos  ■jjy)  for  n  *  0,1,..., 25,  and  the  last  25  data  points 

by  |(1  -  cos  for  n  =  N-25,  N-24, . . .  ,N-1. 

Because  use  of  the  window  appears  to  strengthen  the  detectability 
of  the  multipath  cepstrum  components,  this  window  has  been  included  in 
the  final  version  of  the  computer  program.  However,  the  effect  of  the 
window  on  the  cepstrum  has  i'ot  yet  been  studied  analytically,  nor  an 
optimum  window  shape  determined. 


4.  THE  PSK  SIGNAL 


The  accuracy  with  which  the  multipath  characteristic  of  a  channel 
can  be  estimated  depends  strongly  on  the  type  of  transmitted  signal. 

The  cepstra  of  the  possible  transmitted  signals  have  to  be  sufficiently 
similar  to  each  other,  and  at  the  same  t:me  sufficiently  distinct 
from  the  class  of  multipath  cepstra,  in  order  to  permit  a  reasonably 
accurate  determination  of  the  multipath  characteristic.  Of  course, 
perfect  identification  of  the  multipath  characteristic  can  never  be 
assured  because  of  the  effects  discussed  in  Sec.  3.4. 

4.1  Specification  of  the  PSK  Signal 

The  signals  which  have  been  considered  in  this  study  are  multi¬ 
channel  PSK  signals  as  they  would  be  generated  by  a  16-channel  HF 
data  modem,  *.n  accordance  with  MIL-STD-188C .  This  ;ype  of  signal  has 
a  structure  which  makes  possible  the  extraction  of  useful  multipath 
information  from  the  cepstrum.  It  consists  of  16  differentially 
coherent  phase-shift  keyed  tones  at  frequencies  (933  +  k  •  110)  Hz 
(where  k  *  0,1,..., 15).  Keying  is  simultaneous  in  all  channels, 
at  75  bauds  per  second.  Fig.  12  describes  the  four-phase  encoding 
which  is  specified: 
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odd 


0 

even 


*odd 


0 

even 


0°  (referred  to 

previous  baud) 
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odd 


1 

even 


1 

even 


Fig.  12  Specification  of  Quadra-phase  Modulation 


A  synchronizing  tone  at  825  Hz,  which  is  also  called  for  in  MIL-STD-188C, 
'rtd  an  optional  605  Hz  Doppler  tone,  were  not  incorporated  in  the 
synthesized  signals. 

4.;  1  The  Synthesized  Signal 

Sequences  of  signal  samples,  representing  samples  of  the  PSK  signal 
as  described  above,  were  synthesized.  Independent,  random  modulation 
was  used  on  all  channels.  Such  synthesized  signal  sequences  were  used 
in  all  tests  performed  with  the  computer  pugram  described  in  the  next 
Chapter. 

A  sampling  rate  of  8,175  saraples/sec.  was  used.  This  resulted  in 
exactly  109  samples/baud.  A  signal  sequence  of  512  samples  (or  1024 
samples)  therefore  did  not  consist  of  an  integer  number  of  bauds. 

In  order  to  minimize  the  frequency-domain  aliasing  effect  it  was 
desirable  to  approximate  the  bandlimiting  effect  of  practical  trans¬ 
mission  channels  and  filters.  For  this  reason,  the  synthesized  signaJ. 
was  modified  near  the  baud  transitions  by  a  sine-integral  weighting, 
extending  5  samples  to  either  side  c*  the  baud  boundary.  This  had  the 
same  effect  as  approximately  bandlimiting  the  ideal  lii-chavmel  PSK 
waveform  to  a  frequency  of  4.0875  kHz  (the  Nyquist  frequency)  prior 
to  sampling. 

4 • 3  Cepstrum  of  the  PSK  Signal 

Analytic  derivation  of  the  signal  cepstrum  is  awkward  and  does  not 
lead  to  a  closed  form.  Instead,  the  investigation  has  centered  around 
computer-generated  (discrete)  reel  cepstra  of  the  PSK  signal.  Examples 
of  these  are  shown  in  Figs.  13  and  14.  Only  the  real  cepstra  (inverse 
DFT  of  real  logarithm  of  spectrum)  have  bean  used  because  the  imaginary 
part  of  the  logarithm  tends  to  yield  too  irregular  a  ceustrum.  The 
latter  therefore  has  so  far  not  beer.,  ’tilized  for  multipath  detection. 
This  has  the  disadvantage  that  the  polarity  of  the  multipath  delay 
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2  J 


typical  4096-point  signal  cepstrum 


is  lost.  The  present  processing  procedure  therefore  does  not  identify 
whether  the  weaker  component  in  the  received  signal  precedes  or 
follows  the  stronger  one. 

Some  genera]  features  of  the  real  cepstrum  (figs.  13  and  14)  can 
be  explained.  Large  tluctuations  at  the  first  few  cepstrum  points  are 
the  result  of  the  general  shape  of  the  DFT,  which  is  large  over  the 
range  of  frequencies  where  the  carrier  tones  are  located,  and  falls 
off  to  either  side.  After  taking  the  logarithm  of  the  DFT  magnitude, 
there  results  the  general  form  Bhown  in  Fig.  15.  Additional  distinctive 
features  are  a  peak  corresponding  to  the  channel  spacing,  (at  point  // 75 
ir  Figs.  13  and  14)  and  a  large  ripple  centered  at  the  baud  rate  (point 
//109  in  Fig.  14).  The  remaining  irregular  cepstrum  waveform  is  rela¬ 
tively  small  and  contains  the  detailed  signal  information. 

A  common  component  of  the  various  signal  cepstra  emerges  when  a 
number  of  distinct  signal  cepstra  are  averaged,  as  shown  in  Fig.  16. 

The  degree  of  variability  around  this  common  component  is  illustrated 
by  the  max.  and  min.  plots  in  Fig.  18,  which  are  based  on  a  comparison 
of  nine  distinct  signal  cepstra.  Fig.  17  shows  a  plot  of  the  estimated 
standard  deviation  at  each  cepstrum  point,  based  on  the  same  nine 
signal  cepstra. 


I 


Fig.  15  Sketch  of  general  trend  of  the  discrete  logarithmic  spectrum 
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5.  THE  COMPUTER  PROGRAM 


The  computer  program  for  nultipath  processing  which  has  been 
developed  exists  in  the  form  of  a  FORTRAN  punched-card  deck  of  approx¬ 
imately  550  cards  for  operation  in  the  Honeywell  (formerly  GE)  635/645 
computer  at  RADC.  The  organization  and  use  of  this  program  is  described 
in  this  chapter. 
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Fig.  19  Overall  Functional  Block  Diagram  of  Computer  Processing  System 


5. 1  General  Description 

Fig.  19  is  a  block  diagram  of  the  major  program  functions  as 
described  in  this  Section,  exclusive  of  the  various  plot  and  printout 
features. 

The  program  reads  and  retains  the  contents  of  an  input  data  file. 

It  selects  from  this  data  a  record  of  specified  length.  This  record  con¬ 
stitutes  the  signal  sample  sequence  which  is  processed  during  one  "pass". 
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The  signal  data  undergoes  the  transformation  D  as  described  in  Sec.  2.2, 
with  the  exception  that  only  the  real  part  of  the  complex  logarithm  is 
utilized  for  the  generation  of  the  cepstrum.  The  real  part  of  the  cep- 
strum  is  plotted,  c’T-1  the  values  of  both  the  real  and  imaginary  part  are 
printed  out.  After  this,  there  are  two  possibilities.  If  statistical 
data  is  to  bt  collected  about  the  properties  of  the  cepstrum  (for  use 
in  the  multipath  detection  system; ,  then  a  suitable  routine  can  be  entered 
which  .<  ecords  information  about  the  first  100  points  of  the  cepstrum. 
Thereafter,  the  next  data  sequence  will  be  processed.  After  a  specified 
number  of  passes,  the  mean,  standard  deviation,  min.  and  max.  at  each 
of  the  first  40  cepstrum  points,  will  be  computed  and  printed  ,ut. 

If,  on  the  other  hand,  multipath  processing  is  to  be  carried  out, 
then  the  program  directs  control  to  the  multipath  detection  subroutine. 

It  examines  the  cepstrum  for  the  possible  presence  of  a  multipath  struc¬ 
ture.  If  a  sufficiently  prominent  multipath  effect  is  detected  so  that 
it  can  be  reasonably  clearly  identified,  the  strength  and  time  delay  of 
the  multipath  component  is  estimated  and  this  information  is  printed  out. 
At  this  point  there  exists  the  option  of  terminating  the  "pass",  and 
initiating  the  next  pass.  Otherwise,  the  cepstrum  of  the  calculated 
multipath  component  is  constructed,  and  printed  out,  and  its  FFT  computed. 
The  latter  is  chen  subtracted  from  the  complex  log  of  the  input  signal 
spectrum  so  that,  ideally,  the  contribution  of  the  multipath  is  now 
neutralized.  Subsequ3nt  exponentiation  and  inverse  F.T.  completes  the 
D  ^  transformation.  The  reconstructed  signal  thus  obtained  is  recorded 
into  the  output  data  file.  Control  may  then  be  returned  to  the  beginning 
of  the  program  for  another  "pass",  utilizing  the  next  set  of  data  samples 
in  the  input  data  file. 

Provisions  for  optional  table  print-outs  and  plots  are  also  in¬ 
corporated  at  various  stages  of  processing. 
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5. 2  Data  Files 


As  presently  set  up,  the  program  works  with  several  BCD  data  files. 
These  are  defined  as  permanent  files  emanating  from  the  user's  master 
catalog,  as  follows: 


File  name 

Function 

FI 

input  data  file 

Cl 

output  file 

Ml 

mean  signal  cepstrum 

T1 

threshold  levels 

W1 

weights  used  in  multipath 
detection 

In  order  to  conserve  file  space,  the  files  have  been  written  with  five 
entries  per  record,  using  the  format 

(IX,  E13.6,  IX,  E13.6,  IX,  E13.6,  IX,  E13.6,  IX,  E13.6) 

The  input  data  file  contains  the  signal  sample  to  be  processed. 

The  program  reads  the  entire  input  file,  into  the  array  XIN,  which  is 
presently  set  up  for  1360  sample  values.  A  change  in  this  number  can 
be  effected  by  revising  the  READ  statement  for  FI  and  by  changing  th^ 
dimension  of  XIN  in  the  COMMON  statement  in  the  main  program  and  in  sub¬ 
routines  SPLT,  SPI.OT,  GEN,  and  STAT. 

The  portion  of  the  input  data  which  is  to  be  processed  in  any  one 
"pass"  is  identified  on  the  data  card  by  specifying  the  index  in  the 
array  XIN  which  is  the  beginning  of  the  sample  sequence  to  be  processed. 

The  contents  of  files  Ml  and  T1  are  determined  on  the  basis  of 
measurements  on  the  type  of  signal  to  be  processed.  (These  measurements 
are  performed  in  the  STAT  subroutine.)  The  contents  of  files  Ml  and  T1 
are  read  by  the  DET  subroutine. 

In  file  W1  are  stored  a  set  of  weights,  or  correction  factors, 
which  depend  on  the  data  window  shape  (see  Sec.  5.3.6).  If,  during 


the  detection  process,  a  multipath  cepstrum  spike  is  detected  at 
position  k,  then  its  magnitude  gets  multiplied  by  the  weight  for 
that  positio...  in  order  to  yield  an  improved  estimate  of  the  multi- 
path  gain  a. 

5.3  The  Main  Program  and  Subroutines 

The  detailed  arrangement  of  the  main  program  is  apparent  from  the 
flow  chart  in  Fig.  20.  It  performs  the  operations  described  in  Section 
5.1.  However,  all  the  major  computational  steps  are  carried  out  in  the 
various  subrout  nes.  Communication  between  the  main  program  and  sub¬ 
programs  is  through  COMMON  and  the  subroutine  calling  sequences. 

The  major  program  functions  to  be  executed  are  selected  by 
assigning  the  appropriate  value  to  the  variable  1ST.  This  value  is 
entered  as  data  when  running  the  program.  1ST  has  the  following  effect: 

1ST  *  0  computation  of  cepstrum,  multipath  detection,  and  re¬ 
construction  of  input  signal  with  multipath  removed. 

1ST  =  1  computation  of  cepstrum  and  gathering  of  statistics  on 
the  first  100  cepstrum  points. 

1ST  =  2  computation  of  cepstrum  and  multipath  detection. 

The  program  is  arranged  to  provide  the  following  printer  outputs: 

1)  A  plot  of  the  real  part  of  the  cepstrum,  and  a  listing  of 
the  real  and  imaginary  parts. 

2)  A  listing  of  the  estimated  multipath  cepstrum. 

3)  A  listing  of  the  log-spectrum  of  the  corrected  signal. 

4)  A  listing  of  the  spectrum  of  the  corrected  signal. 

5)  A  listing  of  the  corrected  signal,  under  the  heading 

"BASE  FUNCTION  RESULT". 

In  addition  to  these,  various  choices  of  output  can  be  selected  by 
assigning  the  appropriate  value  to  the  variable  IOUT.  These  choices 
pertain  to 

1)  the  input  data 

2)  the  input  spectrum  (if  plotted,  it  is  the  magnitude  spectrum) 

3)  the  log  spectrum  of  the  input. 
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print  b*se  function  result 


The  value  of  IOUT  is  entered  as  data  when  running  the  program.  IOUT 
has  the  following  effect: 

IOUT  -  i  print  out  tables  of  data  only 

IOUT  =  2  plot  graphs  only 

IOUT  =  3  print  out  tables  and  plot  graphs 

IOUT  =  4  no  output 

IOUT  =  5  data  is  stored  in  file  Cl 

The  main  program  COMMON  statement  specifies  the  following  dim- 
emsioned  arrays: 

X(2, 512, 2) ,  XIN(1360),  XV(512,2),  CP(512),  XM(100,9)  . 

If  data  records  of  more  than  512  samples  are  to  be  processed  in  one 
pass,  the  dimensions  of  some  of  the  arrays  have  to  be  changed  in  all 
the  COMMON  statements  where  those  arrays  occur. 

5.3.1  Subroutine  FFT  -  Fas!  Fourier  Transform 

This  subroutine  is  an  adaptation  of  the  FFT  program  described  in 
[4].  It  performs  both  the  forward  and  inverse  FFT,  as  selected  by 
the  value  of  the  variable  SGN  in  the  following  way: 

SGN  =  -1.  forward  transform 

SGN  =  +1.  inverse  transform 

The  value  of  SGN  is  set  in  the  main  program,  and  in  the  subrouting  DET, 
prior  to  calling  FFT. 

The  subroutine  operates  on  the  data  contained  in  the  array 
X(1,I,J),  and  the  result  is  again  entered  in  X(1,I,J).  The  second 
argument,  I,  identifies  a  particular  (complex)  data  sample,  and  the 
third  argument,  J,  identifies  the  real  and  imaginary  parts: 

J  =  1  real  part 

J  “  2  imaginary  part. 

The  calling  sequence  for  FFT  is 

(N,  NSTAGE,  SGN), 
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where 

N  is  the  number  of  samples  in  the  signal  record  to  be 

processed  (a  power  of  2) 

NSTAGE  is  the  binary  logarithm  of  N 
SGN  as  described  above.. 

The  values  of  N,  NSTAGE  are  read  from  a  data  card  by  the  main  program. 

5.3.2  Subroutine  CLOG  -  Complex  Logarithm 

This  subroutine  computes  the  complex  logarithm  of  the  data  con¬ 
tained  in  the  array  X(1,I,J)  and  returns  the  result  into  the  array 
X(1,I,J).  The  magnitude  of  the  real  part  of  the  logarithm  has  been 
delimited  in  the  negative  direction  r  o  that  it  does  not  exceed  a  value 
of  61nl0  below  the  logarithm  of  VMAX,  where  VMAX  is  computed  in  the  main 
program  as  the  maximum  absolute  value  appearing  in  the  array  X(1,I,J). 
The  imaginary  part  of  the  logarithm  is  computed  in  accordance  with  the 
algorithm  described  in  [5]. 

The  subroutine  also  assigns  values  to  the  variables  XB,  XI,  and  K1 
as  follows: 

XB  argument  of  sample  number  N/2  +  1  (within  CLOG,  this 
variable  is  called  XM) . 

XI  argument  of  sample  number  1  (the  zero-frequency  sample) 
after  XB  has  been  subtracted. 

K1  magnitude  of  linear  phase  component  that  has  been 
subtracted  (see  [5].) 

The  calling  sequence  for  CLOG  is 

(N,  DIV,  XM,  XL,  Kl), 

where 

DIV  is  VMAX  times  10  ^  (assigned  in  main  program) 
and  all  other  variables  have  been  described. 

5.3.3  Subroutine  CALOG  -  Antilogarithm 

This  subroutine  computes  the  antilogarithm  of  the  contents  of 
the  array  X(1,I,J)  and  returns  the  result  into  X(1,I,J).  It  also 
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assigns  the  polarity  of  the  zero-frequency  sample  in  accordance  with 
the  value  that  was  assigned  to  the  variable  XI  in  CLOG. 

The  calling  sequence  for  CALOG  is 

(N,X1) . 

5.3.4  Subroutine  SPLT  -  Plot  Routine  Driver 

This  subroutine  performs  various  chores  in  connection  with  the 
plotting  operation,  and  calls  the  plotting  subroutine  SPLOT. 

When  a  plot  is  called  for  by  the  main  program,  it  transfers 
the  data  to  be  plotted  into  the  array  XV (512, 2)  prior  to  calling  SPLT. 
If  the  plot  is  to  display  absolute  magnitudes  of  complex  numbers,  then 
SPLT  replaces  the  real  parts  contained  in  XV(I,1)  by  the  absolute 
values.  Then  it  prints  the  header  "REAL"  (if  a  magnitude  plot  or  a 
plot  of  the  real  components  is  to  follow)  or  "IMAG"  (if  a  plot  of 
imaginary  components  is  to  follow).  If  a  cepstrum  plot  is  to  be  made, 
it  sets  the  first  and  last  four  samples  of  the  cepstrum  data  equalto 
zero.  (This  is  to  avoid  an  excessive  compression  of  the  ordinate 
scale  of  the  plot,  which  otherwise  results  from  the  very  large  values 
of  the  cepstrum  near  the  endpoints.) 

The  maximum  and  minimum  values  of  the  data  to  be  plotted  are  then 
determined  and  assigned  to  variable,  YMAX  and  YMIN,  respectively.  In 
case  YMAX  s  YMIN,  no  plot  results.  The  number  of  points  to  be  plotted 
is  then  selected  and  assigned  to  NN,  and  the  subroutine  SPLOT  is  called 
to  carry  out  the  actual  plotting  operation. 

The  calling  sequ jnce  for  SPLT  is 

(ISUB,  INN,  IP,  NN), 

where 

ISUB,  INN  and  IP  are  set  in  the  main  program  and  specify  what 
type  of  data  is  to  be  plotted 

NN  this  variable  carries  the  value  of  the  variable  N  in 

the  main  program. 
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5,3.5  Subroutine  SPLOT  -  Plot 


This  subroutine  is  set  up  to  plot  data  on  either  the  high-speed 
line  printer  or  on  the  teletype  terminal.  The  difference  between  the 
two  is  the  line  size,  which  is  set  when  calling  this  subroutine  by 
assigning  either  the  value  101  (line  printe. '  or  61  (teletype)  to 
the  variable  LSZ. 

The  subroutine  selects  an  ordinate  scale  such  that  the  values  of 
YMAX  and  YMIN  fall  within  the  range  of  the  scale.  The  scale  is  selected 
in  such  a  way  that  major  subdivisions  occur  either  at  multiples  of  10, 
of  5,  or  of  2.  The  data  is  then  plotted  to  the  scale  which  has  been 
selected. 

The  calling  sequence  for  SPLOT  is 

(YMIN,  YMAX,  NN,  LSZ) 

where  NN  determines  the  number  of  points  to  be  plotted,  and  is  assigned 
in  SPLT;  and  the  other  variables  are  also  assigned  in  SPLOT,  and  have 
been  described. 

5.3.6  Subroutine  GEN  -  Data  generator 

The  purpose  of  this  subroutine  is  to  generate  a  data  sequence  of 
specified  length  N  for  processing,  from  the  input  data  XIN,  beginning 
at  position  IFL  of  the  input  data.  It  also  applies  a  "Hanning"  window  to 
the  data,  by  weighting  the  first  and  last  25  data  points  by  a  raised 
cosine  curve. 

The  calling  sequence  for  GEN  is 

(N,  IFL) 

where  both  N  and  IFL  are  read  from  a  data  card  at  the  beginning  of  a  pass. 

5.3.7  Subroutine  STAT  -  Statistics  Collector 

This  subroutine  collects  various  statistics  over  several  passes 
through  successive  data  sequences,  during  one  run  of  the  main  program. 
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For  each  of  the  first  100  cepstrum  points,  it  determines 

a)  the  average  value 

b)  the  standard  deviation 

c)  the  maximum  excursion  above  the  average 

d)  the  maximum  excursion  below  the  average 

over  several  passes.  Data  enters  this  subroutine  in  the  array 
XM(100,9).  The  results  arc*  printed  out  in  a  table. 

The  calling  sequence  for  STAT  is 

(ITMQ) , 

where  ITMQ  is  read  from  a  data  card  by  the  main  program  and  determines 
the  number  of  passes,  which  are  to  be  processed.  The  program  is 
presently  dimensioned  to  handle  a  maximum  of  9  passes  during  one  run. 

5.3.8  Subroutine  PET  -  Multipath  Detection 

In  this  subroutine  the  real  cepstrum  of  the  input  signal  sample 
sequence  is  examined  for  the  possible  presence  of  a  multipath  component. 
If  the  presence  of  a  multipath  component  is  determined,  the  magnitude 
and  delay  of  the  multipath  characteristic  are  estimated  in  the  manner 
described  in  Sec.  6.1  and  6.3.  The  detection  procedure  is  presently 
set  up  to  operate  on  the  set  of  cepstrum  samples  beginning  with  number 
2  and  ending  with  number  40.  (The  upper  limit  is  set  by  the  variable 
IRNG.)  The  arrangement  of  the  subprogram  is  described  graphically 
in  the  simplified  flow  chart  in  Fig.  21. 

This  subprogram  produces  the  following  printed  output: 
a)  If  no  multipath  component  is  detected  in  the  specified  range  of 
the  cepstrum,  then  the  output  is  "NO  MULTIPATH  DETECTED".  Also 
printed  are: 

1)  the  index  oi  the  cepstrum  sample  (the  variable  M)  which 
deviated  most  from  tie  recorded  mean  value; 

2)  the  value  C(M),  which  is  the  cepstrum  at  M;  and 

3)  the  amount  by  which  this  valijj  differs  from  the  recorded 
mean  at  M. 


40 


Subtract  mean  values  (file  Ml) 
from  cepstrum  samples 


no 


(no  multipath 
detected) 


Does  the  resultant  data 
exceed  threshold  (file  Tl) 
at  any  sample  position? 


yes 


Locate  position  (M)  and 
lvalue  (Al)  of  most  prominent 
cepstrum  paak 


no 


Is  threshold  exceeded  also' 
at  M+l  or  M-l? 


yes 


Set  estimated  value  of  t 
equal  to  M-l; 
compute  the  estimated  a 
from  value  of  Al  and  the  weight 
for  position  M  (from  file  Wl) 


Compute  estimated  t  and  estimated  a 
on  the  basis  of  relative  magnitude 
of  the  two  strongest  adjacent  cepstrum 
sai  .Aes,  taking  into  account  the 
appropriate  weights  (from  file  Wl) 


' — »-^Print  out  estimated  values  of  t  and 


! 

Construct  multipath  cepstrum 
for  the  estimated  values  of  t  and  a 

1 

Call  subroutine  FFT  to  operate  on  this  cepstrum  | 

1 

fc) 


Fig.  21  Simplified  flow  chart  for  subroutine  DET. 
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b)  If  the  multipath  is  found  to  be  present  with  a  delay  that  is 
essentially  an  integer  multiple  of  a  sampling  interval,  the 
output  is  "MULTIPATH  DETECTED  AT  C(M)".  Also  printed  are: 

1)  the  index  of  the  cepstrura  sample  (the  variable  M)  at  which 
the  largest  multipath  spike  occurs 

2)  same  as  item  (2)  above;  and 

3)  same  as  item  (3)  above 

4)  the  preset  detection  threshold  at  M; 

5)  the  estimated  multipath  delay  x  (the  variable  TU) 

6)  the  estimated  multipath  gain  (the  variable  Al) 

7)  the  synthesized  multipath  cepstrum 

c)  If  multipath  is  found  to  be  present  with  a  delay  that  is  not 
estimated  to  be  an  integer  multiple  of  a  sampling  interval,  the 
output  is  "MULTIPATH  DETECTED  BETWEEN  C(M)  AND  C(M+1)". 

Either  the  +  or  the  -  sign  appears.  M  always  designates  the 
index  at  which  the  greatest  deviation  of  the  cepstrum  from  the 
mean  value  occurs.  If  the  next  sample  exceeds  the  preset 
threshold,  the  +  sign  applies.  If,  instead,  the  preceding 
sample  also  exceeds  the  preset  threshold,  the  -  sign  applies. 

Also  printed  are: 

1)  The  value  of  M 

2)  The  same  information  listed  under  (2),  (3),  (4)  of  (b) ,  for 
both  K  and  M+l 

3)  The  estimated  multipath  gain  if  only  the  cepstrum  at  M  is 
considered  (the  variable  Al) 

4)  The  estimated  multipath  gain  based  on  C(M)  and  C(M+1)  (the 
variable  A2) 

5)  The  estimated  multipath  delay  x  (the  variable  TU) 

6)  The  synthesized  multipath  cepstrum. 

The  calling  sequence  for  DET  is 

(NN.NSTGE) 

where  NN  is  the  variable  which  carries  the  value  of  the  variable  N  in 
the  main  program,  and  NSTGE  has  been  described 
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5.4  Runnira  the  Program 


In  order  to  operate  the  program,  the  signal  data  to  be  processed 
must  first  be  loaded  into  file  FI,  using  the  format  described  in  Sec.  5.2. 
If  the  detection  parameters  which  are  stored  in  files  Ml,  and  T1  are 
applicable  to  the  signal  to  be  processed,  then  the  program  can  be  run 
in  the  detection  mode.  This  requires  only  preparation  of  the  appro¬ 
priate  data  cards  and  appending  them  to  the  program. 

The  first  data  card  bears  only  a  single  integer,  format  (3X,I2). 

This  is  the  value  to  be  assigned  to  ITMQ,  which  is  the  number  of  passes 
to  be  made  in  a  single  execution  of  the  program. 

As  many  data  cards  as  the  specified  number  of  passes  must  then 
follow,  one  data  card  for  each  pass.  Each  of  these  data  cards  bears  the 
values  to  be  assigned  to  N,  NSTGE,  IOUT,  IFL,  and  1ST,  using  the  format 
(3X,I4,3X,I2,3X,I1,3X,I4,1X,I1) .  The  significance  of  N  and  NSTAGE  are 
explained  in  Sec.  5.3.1,  IOUT  and  1ST  in  5.3,  IFL  in  5.3.6. 

If  the  detection  parameters  which  are  stored  in  files  Ml  and  T1 
are  not  applicable  to  the  signal  to  ty  processed,  statistics  on  the 
signal  cepstrum  (without  multipath)  must  first  be  computed.  The  pro¬ 
cedure  for  this  is  the  same  as  described  above.  The  mean  values  for  the 
cepstrum  which  are  computed  in  the  STAT  subroutine  must  be  loaded  into 
file  Ml.  Also,  on  the  basis  of  standard  deviations  and  peak  deviations 
from  the  mean,  as  obtained  in  the  STAT  subroutine,  detection  thresholds 
must  be  established  and  loaded  into  file  Tl. 
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6.  SYSTEM  PERFORMANCE 


6.1  Detectability  of  the  Multipath  Cepstrum 

In  the  present  implementation,  as  described  in  the  previous  Chapter, 
detection  of  the  multipath  cepstrum  (that  is,  determination  of  the 
presence  of  multipath)  consists  of  simple  threshold  detection  of  the 
strongest  cepstrum  value  over  that  portion  of  the  cepstrum  that  is  of 
interest.  This  permits  reliable  detection  of  multipath  over  a  signi¬ 
ficant  range  of  a  and  r.  The  following  discussion  is  based  on  results 
that  have  been  obtained  with  the  synthesized  16-channel  PSK  signal  as 
described  in  Section  4.1,  sampled  at  8175  samples/sec.,  and  processed 
in  512-point  sequences.  With  this  sampling  rate,  a  multipath  delay  of 
less  than  4  ms  will  result  in  the  strongest  multipath  cepstrum  spike  to 
occur  somewhere  within  the  first  34  cepstrum  samples. 

The  detection  characteristic  depends  on  the  threshold  setting.  In 
Sec.  4.3,  the  variability  of  the  signal  cepstrum  was  described  in  terms 
of  measured  standard  deviations,  means,  and  peak  deviations  from  the 
mean.  The  threshold  setting  then  determines  the  probability  of  false 
alarm  (detection  of  multipath  when  none  is  present) .  Only  a  crude 
estimate  of  false-alarm  probability  is  possible,  for  a  given  threshold 
setting,  on  the  basis  of  these  assumptions: 

a)  Successive  samples  of  the  signal  cepstrum  are  characterized  by 
independent  random  variables. 

b)  These  independent  random  variables  are  approximately  Gaussian. 

c)  Additional  variability  due  to  signal  truncation  and  the  effect 
of  the  window  function  is  negligible.  (This  essentially  says 
that  such  additional  variability  is  already  accounted  for  in 
the  model  as  assumed  under  (a)  and  (b).) 

Proceeding  on  the  basis  of  these  assumptions,  suppose  that  the 
threshold  is  set  at 

0  =  u  +  3o  ,  (6-1) 

n  n  —  n 
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where  6  is  the  threshold  at  the  nth  threshold  sample, 
n 

V  is  the  estimated  signal  cepstrum  mean  value  at  the  nth  cepstrum 
n 

sample , 

anis  the  estimated  standard  deviation  of  the  signal  cepstrum  at 
the  nth  cepstrum  sample. 

Let  X  denote  independent  Gaussian  random  variables  with  mean  y  and 
n  n 

standard  deviations  a  ,  with  n  =  a,,...^,  where  d^,d ^  are  the 
limits  of  the  detection  range.  If  the  detection  range  extends  over  32 
cepstrum  points,  then  the  probability  of  false  alarm  becomes 

p£  -  1  -  !P(|Xn  -  uj  <  3o/2 
=  1  -  . 997 33 2  =  0.08. 

This  is  actually  a  pessimistic  estimate,  since  (a)  some  correlation 
between  cepstrum  samples  surely  exists,  and  (b)  available  data 
indicates  that  peak  variations  in  the  signal  cepstrum  tend  to  be  fairly 
restricted.  Of  course,  in  order  to  imr^ve  the  detectability  of  weak 
multipath,  the  threshold  has  to  be  lowered,  which  also  increases  the 
probability  of  false  alarm. 

For  arbitrary  threshold  settings  0n>  the  probability  of  missing  a 
multipath  spike  at  position  n  (assuming  t  to  be  a  multiple  of  the  time 
between  samples)  with  magnitude  a  =  9^  is 

Pm  =  1/2* 
m 

However,  as  discussed  in  Sec.  3.3,  if  t  is  not  an  exact  multiple  of  the 
sampling  interval,  then  the  largest  multipath  cepstrum  component  is 
reduced  in  magnitude,  for  the  same  >.  This  reduc^.on  is  given  by 
sine  v,  where  v  denotes  the  difference  between  x  and  the  nearest  i’Ueger 
multiple  of  the  time  between  samples.  The  mean  height  of  the  largest 
multipath  cepstrum  component,  averaged  over  x,  is  therefore 

f1/2 

I  sine  v  dv  =  -^S i (tt / 2 )  =  0.87. 

-1/2 


Again  using  the  Gaussian  assumption,  the  probability  of  missing  a 
multipath  component  with  a  11  is  therefore  greater  than  1/2.  If 
the  thresholds  zee  set  according  to(6-l),  the  the  multipath  cepstrum 
corresponding  to  a  =  3o  reaches  on  the  average  only  .87  x  3o^  =  2.610^, 
or  .  39o^  below  threshold,  so  that  the  detection  probability  in  the  case 
a  =  30^  is  actually  only  about  0.35. 

It  is  to  be  noted  that  a  "false  alarm"  ana  a  "mis£  can  occur 
simultaneously.  While  a  weak  multipath  spike  may  be  present  that  does 
not  reach  thieshoid,  the  signal  cepstrum  might  exceed  threshold  at  a 
different  point. 

A  reasonable  threshold  function,  based  on  the  sigral  data  in 
Chapter  4,  results  in  a  detection  region  about  as  shown  in  Fig.  22, 
where  the  threshold  curve  has  been  smoothed,  for  easier  viewing.  Also 
entered  in  this  Figure  are  several  contours  of  detection  probability, 
which  have  been  estimated  on  the  basis  of  the  above  assumptions.  The 
rise  of  these  contours  for  large  t  is  due  to  the  degradation  in  the 
multipath  cepstrum  as  a  result  of  truncat  ■>  of  the  input  signal 
(see  Sec.  3.4) . 

There  is  another  effect  which  does  not  appear  in  Fig.  22.  Because 
the  threshold  increases  sharply  for  ver./  small  x,  it  is  possible 
that  the  first  cepstrum  spike  of  a  multipath  characteristic  having 
moderate  a  and  small  x  goes  undetected,  whereas  the  second  spike  is 
detected  since  it  falls  into  a  lower  threshold  region.  Of  course, 
no  useful  multipath  correction  can  be  accomplished  in  that  case,  so 
that  it  can  be  classified  as  the  simultaneous  occurrence  of  a  "miss" 
and  a  "false  alarm". 

6 . 2  Effect  of  Error  in  Estimating  the  Magnitude  of  the  Multipach 

Component 

The  cepstrum  of  a  single  multipath  component  with  delay  t  (r 

rT 

assumed  a  multiple,  of  the  sampling  interval,  x  =  — -  and  relative  gain 


multiples  of  the  sample  spacing 


k  Residual  Cepstrum  Pattern  after  imperfect  cancellation. 


a  is  as  shown  in  Fig.  23.  If  the  estimated  multipath  gain  is  a  f  a 
then  subtraction  of  the  synthesized  multipath  cepstrum  based  on  gain 
o  from  the  actual  multipath  cepstrum  leaves  a  residual  cepstrum  pattern 
which  appears  in  one  of  four  ways  shown  in  Fig.  24 

None  of  these  patterns  represent  a  single  residual  multipath 
component,  but  a  sequence  of  multipath  contributions  with  delays 
x,  2t,  3t,  ...  .  The  same  analysis  applies  to  all  four  cases,  and  is 
conveniently  carried  out  in  terms  of  z-transforms.  Let 


c^  r  denote  the  actual  multipath  cepstrum  at  '.elay  r  samples, 
magnitude  a; 

c^  ^  denote  the  estimated  multipath  cepstrum  ac  delay  r  samples, 


magnitude  a* 
i  of  c  gives 
the  DFT  gives,  in  the  discrete  complex  logarithm  domain: 


Subtraction  of  c  gives  the  residual  cepstrum  c  -  c  Applying 

<?  *r  a  ,r  Q!,r 


cln  (1  +  az  r)  -  cln  (1  +  az  r)  =  cln  ------- 

1  +  az 


Exponentiation  gives 
-r 


1  +  az 


1  +  az 


-r 


1  +  (a  -  a)  z  r  -  a  (a  -  a)  z  a^(a  -  a)z 


(6-2) 


.  (6-3) 


This  series  is  the  z-transform  of  a  transfer  function  consisting  of  a 
sequence  of  impulses,  corresponding  to  multiple  paths,  or  a  tapped 
delay  line  (Fig.  25). 
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As  rough  measure  of  performance,  a  signal- to-interferer.ce  ratio  can 
be  computed  under  the  assumption  that  the  superimposed  delay  components 
are  uncorrelated.  First,  at  the  receiver  input,  a  single  multipath 
component  with  relative  strength  o  gives  an  r.m.s.  signal-to-interf erence 
ratio 


R 


i 


(6-4) 


where 


Vg  =  r.m.s.  signal  voltage 

V  *  r.m.s.  multipath  interference  voltage, 
ra 

After  the  cepstrum  processing  described  above,  there  results  a  super¬ 
position  of  many  residual  multipath  components,  resulting  in  a  m.s. 
value  of 

y.2 .  .2  .2.4.  ,  ..2  (a  -  «)2 

3  -  '  3  1  -  «2  (6-5) 


The  r.m.s.  output  signal-to-interf erence  ratio  is  therefore 

„  _  Vi  -  g2 


1“  *  “I 


(6-6) 


The  improvement  in  signai-to-interf erence  ratio  as  a  result  of 
cepstrum  processing  is  then 


!a  ,£IZ. . 

Ri  |l  -  of/a i 


(6-7) 


In  Figure  26,  contours  of  this  function  are  plotted  vs.  a  and  a. 

R0 

In  the  region  where  —  <  1,  an  actual  degradation  results. 

i 

Howevei ,  the  output  signal-to-interf erence  ratio  alone  may  be  of 

interest.  Thus,  it  may  be  asked,  to  what  portion  of  Fig.  26  must 

operation  be  confined  in  order  that  the  signal-to-interf erence  ratio 

after  processing  does  not  exceed  a  specified  value.  This  question  is 

answered  in  Figure  27.  There,  contours  of  R  =10  (20  db  signal-to- 

o 

interference  ratio)  and  Rq  =  5  (14  db  ratio)are  plotted. 
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Fig.  27  Output  Signal-to-Interference  Ratio  Related  to  Error 
in  the  Estimated  Multipath  Amplitude 


The  effect  of  compensating  for  only  one  multipath  component,  when 
actually  several  multipath  components  exist,  can  oe  analyzed  in  a  similar 
manner.  Suppose  the  received  signal  is  of  the  form  (2-9).  For  sim¬ 
plicity,  assume  t 


riT 


■  j  N  ',  i  =  1,2.  Without  c^rsidering  the  multipath 
cepstrum,  it  is  easily  seen  that  in  this  case  v 6—2)  becomes 


~ri  ~r9  -r^ 

cln(l  +  ctz  +  Bz  )  -  cln(l  +  az  ) 


-r. 


cln 


1  +  az 


±  1?. 


(6-8) 


1  +  az 


Exponentiation  gives 
-r 


1  -r2  rl 

1  +  az  +  Bz  1  +  az 


-r. 


-r. 


1  +  oz 


"rl  _rl 

1  +  az  1  +  az 


(6-9) 


The  first  term  on  the  right  is  the  same  as  (6-3).  The  additional  term 
expands  to 


~r„  -(r  +r0)  2  -(2r  +r  ) 

6z  -  paz  ‘  +  a  8z  - _  (6-10) 


This  shows  that  besides  leaving  the  second  multipath  component  un- 

—  r- 

L  2 

touched  (the  term  8z  ) ,  an  additional  sequence  of  delay  components 
is  introduced  at  delays  of  r^+r2>  ^ri+r2’  etc* 

6. 3  Extraction  of  Multipath  Information 

Assuming  that  a  multipath  component  has  been  detected,  the  error  in 
extraction  of  the  multipath  information  can  then  be  discussed. 

It  has  already  been  mentioned  that  the  present  system  does  not  yield 
polarity  of  t.  This  point  is  considered  further  in  Sec.  7.1. 

Based  on  Equ.  (3-22),  the  present  method  of  estimating  x  consists  of 
fitting  a  sinc-function  to  two  adjacent  cepstrum  samples,  if  both  these 
samples  exceed  the  threshold.  The  center  of  the  sinc-function  then  gives 
the  estimated  value  of  i.  If  the  two  cepstrum  samples  which  straddle 
the  true  value  t  both  fail  to  reach  threshold,  no  multipath  is  detected. 
If  only  one  of  them  exceeds  threshold,  then  the  estimated  multipath 
delay  is  taken  as  the  position  of  that  sample. 

The  error  in  x,  the  estimated  multipath  delay  (aside  from  polarity) 
can  therefore  be  assumed  to  be  restricted  to  within  one  sampling  interval 
of  the  true  value,  giving  |x-x|  <  .122  ms,  and  the  accuracy  improves  with 
increasing  a.  A  precise  statement  is  not  easy  to  formuJate  because  the 
accuracy  of  the  estimate  x  depends  also  on  the  relation  between  x  and  the 
sampling  intervals.  If  x  is  an  exact  multiple  of  the  sample  spacing, 
then  by  definition  of  "detection"  of  the  multipath  cepstrum  spike,  r  =  x. 
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Extraction  of  the  multipath  gain  e  can  be  characterized  graphically 
in  the  manner  shown  in  Fig.  28,  which  applies  for  a  range  of  t  for  which 
the  threshold  is  set  for  approx.  0.1  (see  Fig.  22).  The  probability 
contours  in  Fig.  28  are  based  again  on  the  data  of  Ch.  4,  and  the 
simplifying  assumption  that  the  signal  cepstrum  values  are  governed  by 
Gaussian  r.v. ' s. 

Comparison  of  Fig.  28  with  Figs.  26  and  27  yields  the  probability 
with  which  a  specified  reduction  in  interference  or  a  specified  output 
signal-to- interference  ratio  can  be  achieved. 


Fig.  28  Extraction  of  the  Multipath  Magnitude,  for  r  ■»  2  ms. 


6.4  Reconstruction  of  Transmitted  Signal  Component 


In  the  absence  of  multipath,  straight-through  processing  of  the 
input  signal  via  the  D  and  D  *  blocks  of  Fig.  19  leads  to  an  output 
which  contains  negligible  error.  This  has  served  as  a  check  on  the  sub¬ 
routines  FFT,  CLOG,  and  CALOG,  and  demonstrates  that  there  is  no  signi¬ 
ficant  build-up  of  round-off  errors.  However,  when  using  the  program 
in  the  multipath  detection  mode,  the  reconstructed  signal  appearing  at 
the  output  cannot  be  expected  to  be  "pure",  due  to  errors  in  the 
estimates  of  the  multipath  parameters  (see  Sec.  6.2).  But  even  if  the 
estimate  were  perfect,  serious  distortions  will  arise  near  both  ends  of 
the  output  record  due  to  (a)  the  end  effect  described  in  Sec.  3.6  (see 
Fig.  11)  (b)  the  data  window  (see  Sec.  3.7). 

Therefore,  when  utilizing  the  reconstructed  signal  delivered  by  the 
computer  program  of  Ch.  5,  it  is  important  that  allowance  is  made  for 
these  end  effects.  The  portion  of  the  output  record  which  coincides 
with  the  data  window,  as  well  as  an  additional  guard  space  approximately 
equal  to  the  largest  t  to  be  expected,  had  best  be  discarded.  This 
amounts  to  about  50  points  at  each  end  of  a  512-point  output  sequence. 

A  correction  could  be  applied  at  the  output  to  reduce  the  dis¬ 
tortion  due  to  the  data  window.  This  had  not  been  tried  because  the 
other  effect  (a)  would  scill  remain. 
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7.  DISCUSSION 


7.1  Possible  Improvements  in  the  Present  System 

The  data  presented  in  Chapter  6  shows  thst  with  the  present 
system,  operating  on  a  synthesized  signal,  a  significant  reduction  in 
multipath  interference  can  be  achieved  with  good  reliability,  and  over 
a  wide  range  of  multipath  magnitude  and  delay.  This  statement  assumes 
that  only  a  single  mul.tipath  component  exists,  and  that  the  sign  of  r  is 
known.  The  first  of  these  two  restrictions  can  be  overcome  by  a  more 
sophisticated  detection  routine  which  can  be  implemented  as  a  straight¬ 
forward  extension  of  the  existing  detection  system.  However,  the  same 
type  of  random  error  in  the  estimated  multipath  amplitude,  which  was 
described  in  Chapter  6,  appears  then  independently  on  each  detected 
multipath  component,  so  that  the  probability  of  achieving  a  specified 
reduction  in  the  multipath  interference  is  .  jnd  to  be  less  with  several 
multipath  terms  than  with  only  a  single  multipath  component  (Figs. 

26  to  28). 

Minor  improvements  can  he  achieved  by  optimizing  various  system 
parameters,  such  as  sampling  rate,  length  of  data  sequence  processed  in 
one  pass,  shape  of  the  data  window.  Also,  various  improvements  in  the 
detection  procedure  are  possible.  For  Instance,  the  two  largest 
cepstrum  spikes  (instead  of  only  the  largest  one)  could  be  processed 
in  order  to  achieve  an  improved  estimate  of  the  multipath  amplitude  and 
delay.  However,  such  additional  refinements  (except  means  of  detecting 
the  polarity  of  t)  will  probably  contribute  little  to  clarifying  the 
feasibility  of  the  general  approach  to  practical  multipath  signal 
processing. 

Computation  of  the  phase  spectrum  from  (dx(n)},  using  the  al¬ 
gorithm  described  in  [5],  does  not  yield  a  unique  result.  It  has  there¬ 
fore  not  been  possible  to  obtain  a  meaningful  average  contribution  of 
the  signal  cepstrum,  which  might  then  be  subtracted  from  each  particular 
realization  of  the  phase  spectrum  in  order  to  improve  detectability  of 
the  multipath  contribution. 
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A  different  arrangement  for  computing  the  phase  spectrum  is  pro¬ 
posed  which  yields  a  unique  result.  This  can  be  understood  by  consid¬ 
ering  the  derivative  of  cln  7  (f)  =  ln|  7  (f )  |  +  jt/>  (f ) ,  where  ij;  (f) 


is  the  phase  spectrum. 


j  7(f) 

—  In  7  (f)  -  cy  (f) 


7(f) 

7(f) 


Taking  the  imaginary  part  gives 

„  <(f) 

^x(f)  = 


+  j*x(f) 


(7-1) 


(7-2) 


» 

where  7^(f)  is  the  F.T.  of  (-j2irt  x(t)).  Application  of  the  inverse  F.T. 
to  (7-2)  gives  (in  the  t'  or  cepstrun  domain) 


=  -j2Trt'7"1,  (f) 


(7-3) 


On  the  right  side  of  (7-3)  appears  7  r  ( t * )  ,  which  is  simply  the  odd 
part  of  (?x(t').  This  shows  that  the  even  and  odd  parts  of  the  complex 
cepstrum  can  be  computed  separately.  The  even  part  is  due  to  the  real 
logarithm  of  the  signal  spectrum  and  can  be  computed  in  the  standard 
manner.  The  odd  part  of  the  complex  cepstrum  is  the  contribution  of 
the  imaginary  part  of  the  complex  log  of  the  signal  spectrum.  According 
to  the  above,  it  can  be  obtained  by  taking  the  F.T.  of  —  j  2ti t  x(t) , 
dividing  by  7^(f)»  taking  the  inverse  transform  of  the  imaginary  part, 
and  dividing  by  -j2irt'. 


At  first  it  might  appear  as  if  a  computational  scheme  based  on 
the  above  derivation  introduces  excessive  processing  complexity 
because  two  different  functions  have  to  undergo  Fourier  transformation, 
and  then  inverse  transformation.  This  is  not  so,  however.  Because 
the  two  functions  to  be  transformed  are  real  and  imaginary,  respectively, 
they  can  be  processed  simultaneously  by  the  FFT  algorithm  and  their 
respective  contributions  separated  at  the  output.  Fig.  29  shows  how  the 
complete  D  block  (Fig.  19)  might  be  arranged  to  accomplish  this. 


Of  course,  it  is  not  clear  whether  this  method  of  treating  the 
phase  spectrum  will  be  of  utility  with  the  PSK  signal  of  Chapter  4. 
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Fig.  29  Alternate  Scheme  for  computing  the  complex  cepstrum 

It  will  have  to  be  determined  whether  it  will  result  in  sufficiently 
small  variability  in  the  odd  part  of  the  cepstrum  so  that  useful 
multipath  information  can  be  extracted. 

When  the  communication  signals  are  such  that  reliable  extraction 
of  the  polarity  of  r  is  not  possible,  it  is  possible  to  implement  a 
multipath  compensator  in  the  following  way.  Two  outputs  are  provided, 
corresponding  to  +|t|  and  -|tJ„  The  incorrect  output  can  then  be 
rejected  either  on  the  basis  of  human  inspection  (garbled  messages),  or 
automatically  on  the  basis  of  the  quality  of  a  received  test  waveform. 

7 . 2  Additional  Measurements 

Tests  involving  recorded  signals  from  an  actual  HF  modem,  direct 
and  as  received  over  an  HF  link,  should  turn  out  to  be  informative.  The 
full  scope  of  such  a  test  would  involve  sampling  the  received  signal 
and  recording  the  binary  data  on  taoe.  Portions  of  the  tape  record 
then  have  to  be  converted  into  BCD  format  and  stored  in  the  input  file 
described  in  Sec.  5.2.  ihis  data  can  then  be  processed  and  the  output 
stored  in  the  output  file.  From  there,  conversion  back  to  binary  form 
and  reassembly  into  a  continuous  tape  record  is  possible.  The  re¬ 
constructed  signal  can  then  be  converted  back  into  analog  form  for 
demodulation  in  an  actual  HF  modem. 


Such  a  test  would  incorporate  the  influences  of  many  effects  which 
have  so  far  been  neglected,  such  as: 

Sampling  jitter,  and  distortions  due  to  finite  duration  sampling; 
the  exact  shape  of  the  PSK  waveform  as  generated  in  data  modems; 
the  actual  demodulation  performance  of  a  data  modem; 

the  effects  of  background  noise  (atmospheric,  thermal,  and  quantization 
ncise) ; 

possible  frequency  selectivity  of  the  muJti^ath; 
time  variability  of  the  multipath, 

A  simpler  test  configuration,  and  probably  more  realistic  in  terms  of 
ultima ce  implementation,  would  use  the  computer  to  operate  on  the  3ame  data 
as  described  above,  but  only  for  the  purpose  of  extracting  the  multipath 
information.  The  multipath  data  can  then  be  collected  in  a  file,  and  can 
later  be  used  to  control  the  parameters  of  a  digital  multipath  filter 
for  cancelling  multipath  in  the  test  data.  For  cancellation  of  a  single 
multipath  component,  such  a  filter  would  be  of  the  form  shown  ^n  Fig.  30 
(cf.  [6],  and  other  references  cited 
there).  A  disadvantage  with  this 
digital  multipath  concellation  filter, 
however,  is  the  fact  that  cancell¬ 
ation  can  occur  only  at  delays  which  are 
exact  multiples  of  the  time  between 
samples.  In  other  words,  a  degrad¬ 
ation  arises  due  to  the  quantiza¬ 
tion  of  i.  This  degradation  can 
be  reduced  by  using  a  higher 
sampling  rate  for  multipath  cancellation 

than  fcr  multipath  estimation.  Fig.  30  Recursive  Filter  for 

Multipath  Cancellation 

7 . 3  Possible  Simplifications  of  the  Computer  Program 

The  computer  program  described  in  Sec.  5  was  developed  primarily  as  a 
ressoich  tool,  and  is  not  arrai ged  in  the  simplest  possible  manner.  The 
various  provisions  for  output  printouts  and  plots  add  considerably  to  the 


out 
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compile  and  execution  time,  as  well  as  memory  requirement,  and  can  be 
deleted  in  production- type  runs.  But  there  are  also  simplifications 
possible  in  the  various  signal  processing  steps. 

One  simplification  would  result  in  the  elimination  of  one  FFT.  In 
the  present  system,  the  FFT  subroutine  is  called  four  times  when  multipath 
cancellation  is  performed  (see  Fig.  19).  Instead  of  generating  the  cep- 
strum  of  the  estimatec  multipath  and  then  using  FFT,  it  is  possible  to 
generate  the  exponential  of  the  FFT  of  the  multipath  cepstrum  directly 
(see  Fig.  31a).  This  means  division,  instead  of  subtraction,  to  achieve 
multipath  cancellation,  but  eliminates  the  need  for  the  CALOG  subroutine 
altogether,  as  well  as  one  use  of  TFT.  Furthermore,  the  input  signal  is 
now  stored  and  brought  forward  to  the  cancellation  point  in  the  spectrum 
domain,  instead  of  the  log  spectrum  domain.  Since  the  imaginary  part  of 
the  logarithm  is  not  used  in  multipath  detection  at  present,  this  means 
that  the  complex  logarithm  routine  can  be  simplified  to  carry  ouL  only 
the  log-magnitude  operation. 

Further  simplification  is  possible  if  the  program  is  to  be  used  only 
to  generate  multipaLh  data.  This  results  in  the  block  diagram  of  Fig.  31  b. 

Generate 
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Fig.  31  Simplified  Versions  of  Computer  Processor 
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In  this  Figure  appears  a  block  marked  "simplified  FFT  This  is  meant 
to  suggest  the  possibility  of  a  modified  version  of  the  FFT  for  computing 
only  the  first  40  (or  so)  points  of  the  cepstrum,  which  is  the  only  part 
of  the  cepstrum  being  used  for  multipath  detection  and  extraction.  In  the 
considerably  shortened  overall  processing  of  Fig.  31b,  such  a  modified  FFT 
may  lead  to  a  noticeable  reduction  in  processing  time. 

7.4  Potential  Practical  Application  on  HF  Links 

The  potential  utility  of  multipath  processing  using  the  cepstrum 
depends  on  a  variety  of  factors.  Assuming  that  tests  of  the  type  suggested 
in  Sec.  7.2  turn  out  favorably  and  a  suitable  means  for  identifying  the 
sign  of  t  is  available,  then  there  are  the  following  to  be  considered: 

a)  Expected  rate  of  variation  in  a,t 

b)  How  many  multipath  components  are  likely  to  be  present,  with  what 
distribution  of  strengths 

c)  How  much  improvement  in  the  signal-to-interference  ratio  is 
required  so  that  the  data  modem  produces  an  acceptable  error  rate 

d)  Signal- to-noise  ratio  at  the  receiver. 

These  parameters  determine  the  accuracy  with  which  the  multipath  parameters 
must  be  determined,  at  what  rate  they  have  to  be  measured,  and  whether 
averaging  over  several  measurements  can  be  used. 


Another  important  question  which  has  not  yet  been  considered  is  the 
precision  required  in  the  FFT  and  other  processing  steps,  for  satisfactory 
performance.  A  preliminary  examination  of  FFT  precision,  however,  can  be 
based  on  available  literature.  Thus  Weinstein  [7]  derived  the  following 
estimate  of  the  coefficient  quantization  effect  at  the  output  of  a  radix- 
2  FFT: 


E  -7c 

=  <W6)2  ^ 


(7-4) 


where 


error  variance 
output  signal  variance 


v  =  log2N 

q  =  number  of  binary  digits  in  the  mantissa. 
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An  expression  for  the  fixed-point  roundoff  error  is  also  obtained  in  [7], 
namely , 


where  =  input  signal  variance. 

Suppose  the  combined  effects  of  coefficient  quantization  and  roundoff 

are  to  result  in  an  output  noise-to-signal  ratio  of  no  more  than  -40  db, 

2 

for  N  =  512.  The  input  signal  variance  may  be  set  at  ox  =  1/16  in 

order  to  avoid  overflow,  since  fixed-point  arithmetic  is  assumed.  If 

equa.1  contributions  to  the  output  noise  from  coefficient  quantization  and 

2  2 

from  roundoff  is  assumed,  then  o  /a  =  -43  db  is  substituted  in  both 

Y 

(7-4)  and  (7-5) .  This  results  in  a  requirement  of  9  binary  digits 
(including  sign)  for  coefficient  storage,  and  16  binary  digits  (including 
sign)  for  arithmetic  units.  This  is  nominal  for  existing  FFT  processors 
[81. 

Of  course,  once  a  specific  FFT  processor  design  is  under  consideration, 
a  dpcailed  error  analysis,  similar  to  [9],  can  be  carried  out. 

Combining  some  of  the  suggestions  in  Sec.  7.2  and  7.3,  a  practical 
multipath  canceller  might  take  the  form  of  specialized  hardware  for  estim¬ 
ating  multipath  data  (Fig.  31  b) ,  the  output  of  which  is  used  to  period¬ 
ically  update  the  parameters  of  a  recursive  canceller  (Fig.  30).  This 
canceller  operates  on  the  received  signal  in  real  time,  or  with  a  fixed 
delay  to  compensate  for  the  processing  time.  In  such  an  arrangement, 
processing  speed  can  be  sacrificed,  as  long  as  fne  multipath  structure 
does  not  change  too  rapidly. 
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